Critical scaling near jamming transition for frictional granular particles 
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The critical rheology of sheared frictional granular materials near jamming transition is numer- 
ically investigated. It is confirmed that there exist a true critical density which characterizes the 
onset of the yield stress, and two fictitious critical densities which characterize the scaling laws of 
rheological properties. We find the existence of a hysteresis loop between two of the critical densities 
for each friction coefficient. It is noteworthy that the critical scaling law for frictionless jamming 
transition seems to be still valid even for frictional jamming despite using fictitious critical density 
values. 
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I. INTRODUCTION 

Disordered materials such as granular materials [![, 
colloidal suspensions Q, emulsions, and foams Q have 
rigidity above a critical value of density and do not have 
any rigidity below that value. This sudden change in 
rigidity is known as a jamming transition, which could 
be a key concept for explaining the behavior exhibited by 
disordered materials at zero temperature Q . The coordi- 
nation number changes discontinuously at the jamming 
transition point (point J) for static frictionless spheres 
[^,@. Critical scaling laws, similar to those in the contin- 
uous phase transition, can be used to determine rheolog- 
ical properties of granular particles 0-[I3|- The critical 
exponents of these scaling laws remain a topic of cur- 
rent investigations 04H1- For frictionless grains, where 
inertial effects in the e quat ions of motion are included, 
the present authors pL2Hl4l j used a mean-field theory to 
derive a viscosity rj below point J that satisfies the ex- 
pression ((f> j — 0)~ 4 , as the packing fraction <ft approaches 
the critical jamming fraction </>j from below. We believe 
this result to be asymptotically valid in the hard core and 
elastic limits 

On the other hand, recent studies on the jamming tran- 
sition of frictional grains revealed that elastic moduli, 
coordination number, and density of state are strongly 
affected by the introduction of friction [l5l - f2lj |. How- 
ever, we still do not know the details of the quantitative 
change induced in the rheological properties of the jam- 
ming transition by the introduction of friction between 
grains. Moreover, in spite of many studies on the scal- 
ing law of the jamming transition for sheared friction- 
less particles, there are few such studies on the jamming 
transition of frictional particles. Thus, we need to clar- 
ify whether scaling laws in the vicinity of the jamming 
transition of frictional grains can be used. 

In this paper, we have numerically investigated the 
properties of sheared frictional granular particles in the 
vicinity of the jamming transition. In the next section, 
the details of our numerical results will be presented. In 
Sec. Ill A[ we will explain our setup and models. We will 



demonstrate the existence of hysteresis loops for pressure 
and shear stress in Sec. Ill Bl The three critical area frac- 
tions for sheared granular materials will be estimated in 
Sec. Ill CI where one of them is the true critical fraction 
for the jamming transition and the others are fictitious 
critical fractions. In Sec. Ill Dl we will demonstrate that 
the scaling relations for frictionless particles [l2|, [H[ can 
be used even for frictional systems by using the fictitious 
critical fractions. Finally, we will discuss and conclude 
our results in Sec. Mil 



II. NUMERICAL RESULT 

A. Setup of our simulation 

Let us consider a two-dimensional frictional granular 
assembly. The system includes N grains, each having an 
identical mass m. The position, velocity, and angular ve- 
locity of a grain i are respectively denoted by r^, Vi, and 
u)i. Our system consists of grains having the diameters 
0.7cto, 0.8ero, 0.9cto, and cr 0) where there are N/4 for each 
species of grains. 
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The contact force f ^ consists of the normal part / 



and the tangential part ffj as / 
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The 



normal contact force f\ ■ between the grain i and the 



gram j is given 



by /£> = ^e#)9( ff<j - n 3 )n 
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where h\" and are respectively given by h ™ — 

/c (n) (cr.y - rij) - Ti^v^f and n. l3 = r %J /\r l:j \ for the 

normal elastic constant k^ n \ the normal viscous con- 
stant rf n \ the diameter cr, of grain i, = Ti — Vj, 



Oij = (cj + <Jj)/2 and vffi = (Vi — Vj) ■ n 



Here, 



Q(x) is the Heaviside step function defined by 0(x) — 1 
for x > and Q(x) = otherwise. Similarly, the tangen- 
tial contact force between grain i and grain j is given 

by the equation ff) = min(|^* ) |,/x|/^ ) |)sign(/i^ ) )t ij , 
where min(a, b) selects the smaller one between a and 
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b, and hf) is given by 
the tangential unit vector t 



k^uf} - rPM$ with 

ij {~Vij /\ r ij\: x ij /\ r ij\)- 

Here, few and are the elastic and viscous constants 
along the tangential direction. The tangential velocity 
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and the tangential displacement u\f are respectively 



given by v 

isti 
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i dt v., 

stick «J 



(vi - Vj) ■ tij + (<TiUJi + crjiVj)/2 and 
where "stick" on the integral indi- 



(*) 



cates that the integral is performed when the condition 
\hfj\ < n\f\^ \ or another condition v$vf) < is satis- 
fied mil. 

We study the shear stress S and the pressure P, which 
are respectively given by 



S 



P 



1 

V 



N 



EE' 

\ % j>i 



N 



w (EE' 



An) At) 



An) , M) 

J ij T J jj 



(1) 

(2) 



i j>i 



where V is the volume of the system, and (•) rep- 
resents the ensemble average. Here, we ignore the 
kinetic parts of S and P, which are respectively 

given by S K = - \Ef Pi,xPi,y) / '(mV) and P K = 



Ei Pi ■ p, 



/(2mV), because they are significantly 

smaller than the potential parts in Eqs. ([IJ and ([2]) near 
the jamming transition point. 

In this paper, the shear is imposed along the y direc- 
tion and macroscopic displacement only along the x di- 
rection by the following three methods. The first method 
is the SLLOD algorithm under the Lees-Edwards bound- 
ary condition J22J which we call "SL" for later discussion, 
where the time evolution is given by 



dr 



Pi 

-7T = \-7Vi e x, 

at m 



dPi . 



(3) 
(4) 



with the peculiar momentum r p i = m(vi —jyex) and the 
unit vector parallel to the x-direction e x . In this method, 
the shear rate 7 is a control parameter. 

The second method is quasi-static shearing method, 
which we call "QS" [111 [27| ■ In this method, the shear 
strain A7 is applied by an affine transformation of the 
position of the particles. Then, the particles are relaxed 
under the time evolution equations 



dr l 

dt 
dvj 
dt 



Vi 

m 



3=fi 



(5) 
(6) 



until the kinetic energy per particles becomes lower than 
a threshold value P t h- Then, we repeat applying the 
shear and the relaxation process. Here, we chose A7 = 



1(T 6 and P t h = Kr 7 fc(™)cr2, which are small enough not 
to influence our results. This method is expected to cor- 
respond to the low shear limit of the SL method. 

The third method is stress control method, which we 
call "SC". Here, the time evolution equations are given 
by Eqs. (J3]) and under the time evolution of the shear 
rate 7 



dq Sp-S 
~di ~ Q 



(7) 



with the relaxation constant Q = 10 4 fc(")/ero [22j |. Here, 
we choose So = 10~ 7 , which is small enough that the 
results do not depend on the value of So- 

In this paper, we mainly use the SL method for the 
analysis. We also examine the QS and SC methods in 
Figs. [6l and [10] in Sec. II. C to determine the true 
jamming density. 



B. Hysteresis in rheological properties 

First, let us consider rheological properties of sheared 
granular particles. Here, we study the shear stress S and 
the pressure P by using the SL method. We find that 
the shear stress S or the pressure P suddenly changes 
for a finite friction constant fx at a critical value of 
7 near the jamming point. This might be related to 
the existence of a hysteresis loop for frictional granular 
materials pil [28|. To verify the existence of the hys- 
teresis loop, we first vary the shear rate 7 from 70 = 
5.0 x 10~ 4 \/£;(")/m to sequentially decreasing values as 
7 = 70, a _1 7o, a~ 2 7o> ■ • • j a _JVs 7o with the rate of change 
a and the number of the step N s = 2 In 10 / In a, where the 
smaller 7 is 5.0 x 10~ 6 -\/ k^> jm. We call this process the 
decreasing process. Next, we vary the shear rate from an 
initial values of a~ Ns jo to sequentially increasing values 



given by 7 = a N "j ,a N » lA f , ■ ■ ■ ,a 2 7 ,a M'tbTo- 
We call this process the increasing process. Here, we 
vary the shear rate so the data as to be apart in the 
logarithmic plot. For each value of shear rate 7, the sys- 
tem attains a quasi-steady state after the waiting time is 
larger than r w = O.47 -1 , and we have measured the shear 
stress S and the pressure P during the time Tj = 2.07 -1 
after the waiting time is equal to r w . 

As expected, there exists a hysteresis loop in frictional 
jammed systems under this set up. Figure Q] shows the 
scaled peculiar momentum Pi/y/ mT with granular tem- 
perature T = ^2i\Pi\ 2 /(2mN) in both the decreasing 
process (Fig.l (a)) and the increasing process (Fig.l (b)) 
for n = 2.0, cj) = 0.795 , N = 8000, a = lO" 05 and 
7 = 7.9 x 10 -6 -\/fc( n ) jm. In the decreasing process (Fig. 
[IJa)), the motion of grains is localized and the charac- 
teristic length scale is significantly larger than the grain 
size. In the increasing process, as shown in Fig. [TJb), the 
heterogeneity of grain motion is rather suppressed. It 
should be noted that the coordination number Z in Fig. 
[IJa) is larger than 3, whereas Z in Fig. [IJb) is less than 
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2. This difference in the coordination numbers might be 
related to the hysteretic behavior of factional granular 
particles. 

The shear stress S and the pressure P exhibit clear 
hysteresis loops, as shown in Fig. [2j In Fig[2ja), we 
plot the shear stress S in a quasi-steady state of each 
shear rate 7 with /z = 2.0, N = 8000, and a = 10° 1 for 
packing fraction values of <f> = 0.870,0.810,0.795,0.790, 
and 0.780. For highly packed systems such as c/> = 0.870 
or 0.810, 5 or P becomes a constant in a weak shear 
limit (7 — > 0), which implies the existence of the yield 
stress. In contrast, S and P are proportional to j 2 for 
a relatively low packed system at <f> = 0.780. It should 
be noted that S and P are independent of the process 
for both limits, as in the case of frictionless particles 
i, GJ, G3. However, S and P depend on the pro- 
cess for intermediate packed systems at 4> = 0.795 and 
0.790, and thus, hysteresis loops appear in this region. 
For <j> — 0.795, S and P suddenly decrease at around 
7 = 5.0 x 10~ 6 -\//c( ,l )/m in the decreasing process, but 
increases around at 7 = 5.0 x 10 -5 \J M") Jm in the in- 
creasing process as shown in Fig. [2ja). The upper branch 
in the decreasing process is called the solid branch and the 
lower branch in the increasing process, the liquid branch. 

Let us check how the hysteresis loop, as shown in Fig. 
[2] depends on the system size. In Fig. [3J we exam- 
ine the scaled shear rate 7* = jy/k^ jm dependence of 
the scaled shear stress S* = S/ik^cr^ 1 ) for <j> = 0.795, 
fi = 1.8, and a = 10 01 with N = 8000 and 16000. Al- 
though the critical 7* where the two branches appears 
slightly differ between N = 8000 and 16000, the shape of 
the hysteresis loop is almost unchanged. Hence, we con- 
clude that the hysteresis loop still survives even in the 
thermodynamic limit (N — > 00). 

In Fig. [U we compare the hysteresis loop of S for 
a = 10 - 1 with that for a = 10 05 . Owing to the slow 
change of the shear rate, the transition from the solid 
branch to the liquid branch for a = 10 05 takes place at 
a larger shear rate than that for a = 10 . However, the 
existence of the hysteresis loop and the values of S in 
the solid and liquid branches are unchanged even when 
we use a smaller a. Thus, qualitative behavior of the 
hysteresis loop is insensitive to the choice of the change 
rate of 7* . 

Finally, we check the dependence of the hysteresis 
loops on the waiting time r w for each shear rate. In 
Fig. [U we compare the hysteresis loop for r w = O.47 -1 
with that for r w = 4.O7 -1 and r w = 8.07 -1 . Although 
the area of the hysteresis loop for large waiting times is 
smaller than that of short waiting time, the area is con- 
verged for r w is larger than 4.07 -1 . Thus, the hysteresis 
loop survives even when we slowly change the shear rate. 

C. Critical densities and phase diagram 

In this subsection, we demonstrate the existence of 
two critical densities. First one is the transition den- 



sity <j>c{n) & t which the shear stress S in the low shear 
limit 7 — > has a finite value. In order to determine 
4>c(.n)i we introduce the jammed fraction / obtained 
from the simulation using the QS method, where we in- 
troduce / as a fraction of samples where S is larger than 
a threshold value Stb = 10~ 7 £;( n ' /a Q . In Fig. [51 we 
plot the jammed fraction as a function of the area frac- 
tion (f> for /i = 0.0, 0.2 and 2.0. The jammed fraction 
suddenly changes around 4> = 0.843, 0.825 and 0.795 for 
fi = 0.0,0.2 and 2.0, respectively. Then, we can deter- 
mine </>c(m) as the area fraction at / = 0.5. 

The second critical density is obtained from the pres- 
sure in the low shear limit. In Fig. we plot the value 
of the pressure P and the stress S obtained from various 
methods SL, QS and SC as a function of the area frac- 
tion 4> for /i = 0.0 and 2.0. As shown in Fig. [7J P and S 
converge to the value of the QS method in the low shear 
limit of the SL method. For the SC method, if we take 
small enough So, the value of P is almost equal to that 
of the QS method as shown in Fig. \7\ while S becomes 
almost zero. Here, we should note that the pressure P 
and the stress S continuously increase from the transition 
point </>c(/-*) for [i = 0.0, while P and S in the low shear 
limit discontinuously change at 4>c{v) f° r A* — 2-0. This 
discontinuous change should be related to the emergence 
of a hysteresis loop. Therefore, </>c(m) can be regarded as 
a true critical area fraction for the jamming transition. 

On the other hand, we find that a critical scaling is 
satisfied if we introduce another critical point 4>s{lJ-) as 

P(<f>, n) = ILp((f>-<l>s<Ji)), S(<f>, n)/A(ji) = n s (0~0 s ( M )) 

(8) 

for ns(n) < 4>c where the denominator A(fi) in the 
second equation depends only on //. Figure [5] clearly ver- 
ifies the validity of the scaling relation ([5]), though the 
data for c/> < 4 > c(l^) are not involved. Here, we choose 
0s(O) = 4>c(0) for M = 0, while we choose <f>s((J, from 
the collapse of the data onto the universal scaling curve 
for /1 > (Fig. [H]). Our result also suggests that Hp(x) 
and Us(x) are linear functions of x for x > 0. In the 
inset of Fig. [SJ we enlarge the region near (ps((J,) — <j> = 0, 
which indicates that P and S discontinuously change for 
fi > 0.8, while they seems to have continuous changes for 
fx < 0.4. Here, we should note that (ps(n) is a fictitious 
critical fraction, at least for fi > 0.8, but is important in 
characterizing scaling laws for the rheology of the fric- 
tional granular materials in the next subsection. 

In Fig. [9l we show the dependence of the critical frac- 
tions 4>c{n) and </>s(/z) on the friction coefficient p,. Both 
critical fractions decrease as /1 increases. This result is 
similar to the previous ones for three-dimensional fric- 
tional spheres [15l - t2l| . As shown in Fig. [9l the difference 
between <^c(a*) and <f>s(lJ>) is not visible for \l < 0.4, but 
it becomes distinct as 4>c{v) > </ ) s(m) f° r A* > 0-4- This 
discrepancy between ^cCa*) an d ^(m) is related with the 
discontinuous change of P shown in Figs. [7]and[8l Fig- 
ureOalso exhibits the region where the hysteresis appears 
for N — 8000 and a — 10 01 . Here, let us introduce a to 
represent the area of the hysteresis loop shown in Fig. 
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FIG. 1: Snapshot of the scaled peculiar momentum field Pj/v mT for the decreasing process (a) and the increasing process (b) 
for u = 2.0, <j> = 0.795, TV = 8000, a = 10° \ and j = 7.9 x 10~ 6 s/k( n ) Jm. 
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FIG. 2: (a) Scaled shear stress S 1 * = S/(k^ n '<7 Q as a function of scaled shear rate 7* = jy / k("'> Jm for various values of 
packing fraction cf> for jit = 2.0, TV = 8000, and a = 10 01 . 

(b) Scaled pressure P* = P/ (k^- 71 ' Uq 1 ) as a function of 7* for various values of packing fraction 4> for fi = 2.0, TV = 8000, and 
a = 10 01 . 



Ha) 

a = /" " dG {log 10 S+ (7) - log 10 5_ (7)} , (9) 

id 

where 5+ (7) and >5_(7) respectively represent the shear 
stresses in the decreasing and the increasing process for 
G = log 7, Gi = log 10 (70a" JV »), and G 2 = log 10 (70). In 
Fig. |9l we plot a in the [i-<j> plane with a gray scale. The 
hysteresis only appears in a restricted density region for 
each values of fi, and the range of the region becomes 
wider as the friction coefficient /1 increases. In Fig. |H1 we 
also enclose the region of the hysteresis loop by thin solid 
lines, where a is larger than a threshold value ath =0.1. 
The boundary lines are almost identical to 4>c{p) and 
4>s(p) and the hysteretic region lies between 4>s(p) and 



</>c(a*), although a discrepancy from the critical densities 
exist. We suppose that the boundary lines coincide with 
4>c{lA a- n d 4>s{p) if we take a limit a t h — ► 0. However, 
if we use smaller a t h in the present numerical data, we 
could not draw the clear lines due to sampling error. In 
this figure, the hysteresis loop is only verified in the re- 
gion \i > 0.5, though we cannot exclude the possibility 
of the existence of a small hysteresis loop for smaller \i. 
We note that Fig. [§] involves the other fictitious criti- 
cal fraction </>i(/i), which will be introduced in the next 
subsection. 

Let us check the finite size effect of the critical fraction 
4>c{lA- In Fig. I10[ we examine the jammed fraction / 
as a function of for fi = 2.0 with TV = 1000, 2000 and 
4000. Fig. [TU1 indicates the existence of a finite size effect, 



5 



10"' 



10"' 



-4 



10 

S* io" 5 

10" 6 



10" 
10" 1 







II 


f T 








/ - 




f N = 8000 ■ - 




N = 16000 o 



10"' 



10" 



10"' 



10" 



7 




FIG. 3: Scaled shear stress S* = S/(k^a 1 ) as a function FIG. 5: Scaled shear stress S* = <Sy(fc'™Vo) as a function 



of scaled shear rate 7* = -yyfeW jm for <^> = 0.795, jj, = 1.: 
and a = 10 01 with AT = 8000 and 16000. 
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FIG. 4: Scaled shear stress S* = S/(k (n) ao ) as a function of 
scaled shear rate 7* = -y^feM /m for = 0.795, /1 = 1.8 and 
N — 8000 with the shear rate step a — 10 01 characterized 



by the solid squares and a - 
squares. 



10 u 



characterized by the open 



of scaled shear rate 7* = 7 Jm for <f> = 0.795, ^1 = 1.8, 
a = 10° 1 and JV = 8000 with the waiting time r w = 0.4-y" 1 
characterized by the open squares, 4.07 _1 characterized by 
the solid squares and 8.07 _1 characterized by the open circles. 




FIG. 6: Jammed fraction / as a function of <j) f° r N = 4000 
with n = 0.0, 0.2 and 2.0. 



where 4>c{n) lies between 0.794 and 0.798. However, this 
finite size effect does not affect the qualitative behavior 
of the phase diagram in Fig. [5] because the difference 
between <j> = 0.794 and 0.798 is so small that it cannot 
be distinguished in Fig. 



D. Critical scaling laws 

In a frictionless system with the linear spring repulsion, 
the shear stress S 1 , the pressure P and the coordination 
number Z for cf> > <j)j satisfy S oc cf> — <j)j, P oc <f> — 
<j>j, and Z — Z c oc (4> — (jjj) 1 ^ 2 , where Z c = 4 is the 
coordination number for the isostatic state of frictionless 
particles [HI, [l3[ . Remarkably, these scaling laws are still 
valid even for frictional systems if we use <f>s (/x) defined 



in Sec. Ill CI as 



S ~ (f)-cps(lj), P ~ i>-<t>s(fj), Z-Z c (fl) - y/<f>-<f>s(jt)- 

(10) 

Figure [5] supports the scaling relations (|TU1) for the pres- 
sure P and the shear stress S. Here, it should be noted 
that the scaling relation P ~ (<f> — 0,/) 108 , which is re- 
ported in Ref. (HJ, might be better than Eq. (fTQ|) for 
/x = 0.0, but the data for /i > suggests that the linear 
relationship between P and <p — </>s(/u) might be better 
than P ~ (4> - c (/i)) 108 . In Fig. El we plot Z as 
functions of <f> — 4>s(fJ>), which also supports the scaling 
relation (flU)) for Z. The inset in Fig. [11] shows the [i- 
dependence of Z c , which approaches the critical value 
(Z c = 3) for isostatic frictional grains. This behavior 
of Z c (/j,) is similar to the previous result obtained for the 
static granular packing problem f]~5l— f]~TL f]~9 . 21]. It should 
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FIG. 7: (a):The rescaled pressure P* with P* = P/(fc (n) cr l J 1 ) as a function of (j> for AT = 8000 by using the SL method with 
7* = 5 x 10" 5 which we call "SL1", the SL method with 7* = 5 x 10~ 6 which we call "SL1", the QS, and the SC methods for 
fi — 0.0. (b):The rescaled pressure P* with P* = P / (k^ o^ 1 ) as a function of cj> for N = 8000 by using the SL method with 
7* = 5 x 10" 5 which we call "SL1", the SL method with 7* = 5 x 10~ 6 which we call "SL2", the QS, and the SC methods 
for n = 2.0. (c):The rescaled stress S* with 5"* = S/(fc (n) a _1 ) as a function of (j> for N = 8000 by using the SL method with 
5 x 10" 5 which we call "SL1", the SL method with 7* = 5 x 10" 6 which we call "SL2", and the SC methods for fi = 0.0. 
Note that we eliminate the data by SC. (d):The rescaled stress 5" with S* = S / (k^ o^ 1 ) as a function of (f> for N = 8000 by 
using the SL method with 7* = 5 x 10~ 5 which we call "SL1", the SL method with 7* = 5 x 10 -6 which we call "SL2", and 
the SC methods for [i = 2.0. 
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be noted that only the data in the solid branch satisfy the 
scaling relation Eq. ([TU]) . 

For <j) < §J in a frictionless system it is known that 
Bagnold's scaling holds, and S and P satisfy S oc j 2 {(j)j — 
4>)~ A and P oc 7 2 (0j — <fi)~ 4 . Even in the frictional sys- 
tem, S and P satisfies Bagnold's scaling S oc ■y 2 and 
S oc 7 2 in the liquid branch as shown in Fig. [2j but 
any scaling relation using 0c(/-O and <l>s{l*) cannot be 
applied. However, by introducing another fictitious crit- 
ical fraction 4>l{h) shown in Fig. [HI S and P satisfy the 
scaling relations 

S ~ 7 2 {<M/i)-<«- 4 and P~ 7 2 {<M^)-0r 4 (H) 

in the liquid branch. In Fig. we plot S/j 2 and P/j 2 
in the liquid branch, where Eq. is satisfied. However, 



4>l{^) is just a fitting parameter for Eq. (jlip . and could 
not be estimated from an independent protocol. In this 
sense, Eq. (ITT]) might be superficial. 

We should note that the scalings of Z, P and S in the 
jammed phase of the static granular packing problem are 
not affected by the introduction of the friction [l8l . ll9l . l2lj ]. 
Equations (fT0|) and (ITTj) extend the validity of such an 
idea even for a system in the sheared granular systems. 

It should be noted that the range of the excess density 
for the scaling of the liquid branch in Fig. [T^] is at most 
one decade because the shear rate should be smaller to 
obtain the liquid branch near the critical point. However, 
the validity of the scaling in Eq. (jTTj) for the frictionless 
system in the hard core limit has been already verified in 
Ref. 03 in details. 
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■ s /kMjm and AT = 8000 as a function of <A 
on /i i.e., A = 1.0,0.4,0.3,0.3, and 0.3 for /i = 0.0,0.2,0.4,0.8, and 2.0, respectively. Inset shows S* /A near <f>s(n) 
fi = 0.2,0.4,0.8. 
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FIG. 9: The critical fractions <j>c((J,), 4>s(n), and 4>l{h) as a 
function of /i. The amount a of the hysteresis-dependence in 
H — 4> plane is plotted with a gray scale. The region of the 
hysteresis loop, which is defined by the region where a > 0.1, 
is enclosed by thin solid lines. 



III. DISCUSSION AND CONCLUSION 

Let us compare our results with those of the previ- 
ous studies for the static granular packing problem. In 
Refs. (l5l - [2lj |. the ^-dependence of the scaling laws and 
the critical point are studied, where the critical fraction 
decreases as the friction coefficient /1 increases. These re- 
sults are consistent with our results for the solid branch. 
However, the results corresponding to those for the liquid 
branch in Figs. and have not been reported in any 
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FIG. 10: Jammed fraction / as a function of <f> for fi — 2.0 
with N = 1000, 2000 and 4000 obtained from the simulation 
using the QS method with A7 = 10" 6 and E th = 10" 7 fc (n) crg. 



previous studies. Thus, we stress that the scaling laws 
with the aid of the ^-dependence of <?!>l(/x) are our new 
findings. 

Moreover, similar phase diagrams to Fig. |9] are pro- 
posed in Refs. 0, [3l|. However, their results suggest 
that critical densities continuously exist between two crit- 
ical values, while our results indicate that the critical 
densities are at most three. This is the main difference 
between the phase diagram in Refs. 0, HH and ours. 

In Ref. [3(| , it has been already reported that the crit- 
ical fraction splits for frictional granular particles driven 
by a constant force from the boundary wall. There are 
three critical fractions which are denoted by cf>j 1 , cf>j 2 and 
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of — 4>s{n) in the soZid branch. Inset shows Z c (fJ,) as a 
function of the friction coefficient //. 



0j 3 . Here, </>ji is the critical fraction used to characterize 
the appearance of the jamming transition, <f>j2 is that at 
which the viscosity r\ diverges, and (f>j3 is that from where 
the pressure increases. Thus, we expect that these criti- 
cal fractions 4>j2, and <f>j3 respectively correspond to 
<f>c((j), 4>l{^), and 4>s{^) in our system. However, these 
critical fractions satisfy <f>j\ < <j)j2 < <j>j3, in contrast to 
the relation <fis((f) < 4>c{lA < <Al(m) We plan to clarify 
the origin of these differences in the future work. 

In Ref . [lH , the scaling relations are plotted as a func- 
tion of the distance from isostaticity Z — (D + 1) with the 
dimension D of the system, where D + 1 is the coordi- 
nation number of the isostatic state. However, we adopt 
the scaling relation as a function of the distance from the 
critical density <f>— <fis(n) m Fig. [TT] This is because Z is 
not a control parameter and the scaling using Z—(D + 1) 
cannot be used for the liquid branch, where Z tends to 
zero for sufficiently hard grains. 

In previous studies for sheared frictionless materials, 
the different scaling from Eq. (fTU|) are proposed @, HHUI . 
In particular, S ~ (<f> — 4>.i) 3 ^ 2 for the jammed granular 
phase is used in Ref. [l(J, LLLI ■ It should be noted that 
the difference between ours and Hatano's exists only 
in this relationship, but there are several differences be- 
tween Hatano's and the exponents by Tighe et al. 

We also note that there are no similarities between 
Tighe et al. [l0| and Olsson and Title Q- Fortunately, 
our scaling relations S cx <fi — 4>s{^) and P oc <f> — 4>s(n) in 
the jammed region is supported by the simulation pre- 
sented in Figs. [8] and 10. Moreover, the exponent is 
close to the exponent 1.08 obtained by Olsson and Ti- 
tle [Hj], although they claim that 1.08 is an evidence for 
non-linear behavior. It should be noted that the scaling 
relations in Refs. [l2l - fl4j are obtained for granular parti- 
cles with inertia effect, and the scaling relations are held 
in hard-core and elastic limit. Therefore, our scaling ex- 
ponents are not necessary to be identical to those in Ref. 



7, 10] for sheared foams without inertia effect. Indeed, 
the derivation of our exponents is due to the propagation 
of phonon under the isostatic condition (33J. On the other 
hand, HatanofljJ indicated the estimation of the critical 
exponent depends on the choice of the value of the criti- 
cal density <j>j and he suggested P oc (0 — from his 
simulation. However, as shown in Sec. Ill D[ if we chose 
the value of 4>s (/i) which is estimated from the data under 
the three different methods, our scaling exponent seems 
to be better than his [llj. Here, we should note that the 
crossover from P ~ {<f> — <f)j) to P ~ (4> — (f>j) 3 ' 2 might 
be understood from the deviation from the critical point 
as P ~ Z(<j))((j) - <t>j) with Z(<t>) - Z c oc (0 - </>j) 1/2 - We 
will have to resolve the contradiction between ours and 
his in near future. We finally note that the data reported 
in Ref. [lj| also suggests that P oc (<f> — </>j) 3 / 2 , but his 
data is obtained from the Hertzian contact model, where 
the exponent 3/2 is equivalent to ours [13, EH ■ 

Here, we should discuss the macroscopic friction coeffi- 
cient S/P. In Figs. [Tl?land[T4l we plot S/P as a function 
of the density <f> for both the solid and the liquid branch. 
S/ P is almost independent of 4> in the solid branch, while 
the apparent ^-dependence of S/P in the liquid branch 
is observed for the off-critical region <f> — 0l(m) > 10~ 2 , 
which might be related to the discrepancy of S and P 
from the scaling law given by Eq. (fTTj) in Fig. [T^] be- 
cause we assume that S/ P is independent of 4> when we 
derive Eq. (JTTJ) in Ref. [II]. As shown in Fig. HU S/P 
is almost independent of the friction coefficient fi except 
near \i = 0. The /x-dependence of S/P will be discussed 
elsewhere. 

The existence of hysteresis is also known for frictionless 
granular particles under a finite shear stress [3~jj . How- 
ever, this hysteresis differs from that of our frictional sys- 
tem shown in Fig. O In Ref. (34|, the stress is a control 
parameter, and the states with 7 = and 7 / are 
coexistent for a given shear stress. On the other hand, 
the shear rate 7 is a control parameter in the SL method 
used for Fig. [2 and the two states for a given 7 are those 
with large and small shear stress. 

In conclusion, we numerically investigate the sheared 
frictional granular particles, and find the existence of the 
hysteresis loops for different values of the pressure and 
the shear stress, whose relevancy is numerically verified. 
It is confirmed that the critical densities which charac- 
terize the jamming transition are split into three values, 
where one of them is the true critical density for the jam- 
ming and the others are fictitious critical densities. It is 
also verified that the scaling relations (fTU|) and JTTj) for 
frictionless particles can be used for the frictional systems 
by using the fictitious critical densities. 
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